!********************************************************************************************
!       THIS PROGRAM ACCOMPANIES Fried et al (2024)
!
!       "Understanding the Inequality and Welfare Impacts of Carbon Tax Policies"
!
!   This code is used to calibrate the benchmark model
!   If you use this program, I would kindly ask that you cite my paper.
!*********************************************************************************************

PROGRAM hetero_ss
USE functions_mod
USE grid_interpolation
USE lambda2_optimizer_mod
USE retired_valuer_mod
USE working_valuer_mod
USE parameters_mod
USE parameters_energy_mod
USE NEQNF_INT
USE qsort_c_module


IMPLICIT NONE
INTERFACE
    function randu(iseed)
        INTEGER, intent(inout)   ::iseed
    DOUBLE PRECISION   ::randu
end function randu
END INTERFACE
include 'mpif.h'

parameter send_data_tag = 2001, return_data_tag = 2002

INTEGER            :: con_pos,kstart, capital_positive, tax_loop, non_clear,temp,kit,i,s,it1,intloc1
INTEGER            :: ac,out_loop,it,e_it,parm_it,parm_it2,parm_it3,parm_it4,parm_it5,parm_it6,parm_it7,krow
INTEGER            :: sim_shocks(J,simulationsN),sim_types(J,simulationsN)
INTEGER            :: taxes_clear(max_outer_loops+1)
INTEGER            :: shocks_start1,shocks_start2,parm_loop
INTEGER            :: num_procs_run,seed,it2

INTEGER            :: my_id, root_process, ierr, status(MPI_STATUS_SIZE)
INTEGER            :: num_procs, an_id, e_c_l,e_c_h
INTEGER            :: num_procs_run2
INTEGER            :: start_row,end_row,increment,large_inc,increment2,large_inc2
INTEGER            :: bad_combination_inner_nostop,bad_combination_outer,bad_combination_inner_stop

DOUBLE PRECISION   :: chi_l,chi_h,beta_l,beta_h,chi_grid(chiN),beta_grid(betaN),rate_use
DOUBLE PRECISION   :: alevel_grid(aN),amult_grid(amultN)&
    ,kgridl_l,kgridl_h,kgridl_grid(kgridlN)
DOUBLE PRECISION   :: theta_l,theta_h,lambda_l,lambda_h,theta_grid(thetaN),lambda_grid(lambdaN)
DOUBLE PRECISION   :: parms(betaN*chiN*thetaN*lambdaN*1*1*kgridlN,7),xguess(3),solver_out(3)
DOUBLE PRECISION   :: b,pricefactor,energy_share_sim(J,simulationsN),K1_share,N1_share
DOUBLE PRECISION   :: curv,diff_parameter,diff_parameter_v(6),survive(J), dying(J),value,types(typesN)
DOUBLE PRECISION   :: growth,  pop_temp1(J),hours_jnr,sim_earnings(J,simulationsN),consumption_temp_sorted(simulationsN)
DOUBLE PRECISION   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),sim_consumption(J,simulationsN),energy_temp(simulationsN),consumption_temp(simulationsN),consumption_high(J),consumption_low(J)
DOUBLE PRECISION   :: sim_energy(J,simulationsN),energy_high(J),energy_low(J)
DOUBLE PRECISION   :: agg_labor(J,typesN),agg_capital(J+1,typesN),agg_consumption(J,typesN),agg_energy(J,typesN)
DOUBLE PRECISION   :: kprime(typesN,shocksN,kgridN,J-1),labor(typesN,shocksN,kgridN,Jnr)
DOUBLE PRECISION   :: dummy,borrowers,negative_value,negative_count
DOUBLE PRECISION   :: ss(typesN,(shocksN)/2),ss1(typesN,(shocksN)/2)
DOUBLE PRECISION   :: v(typesN,shocksN,kgridN,J+1),v_retired(typesN,shocksN,kgridN,J+1)
DOUBLE PRECISION   :: v_in(typesN,shocksN,kgridN),lambda0_1
DOUBLE PRECISION   :: ctot80(kgridN,1),e80(kgridN,1),c80(kgridN,1), con(kgridN)
DOUBLE PRECISION   :: aggregates(betaN*chiN*thetaN*lambdaN*1*1*kgridlN,48),con_dummy
DOUBLE PRECISION   :: aggregates_ss(betaN*chiN*thetaN*lambdaN*1*1*kgridlN,typesN,(shocksN)/2)
DOUBLE PRECISION   :: aggregates_ss_earnings(betaN*chiN*thetaN*lambdaN*1*1*kgridlN,typesN,(shocksN)/2)

DOUBLE PRECISION   :: energy_miss,N1_ss(typesN)
DOUBLE PRECISION   :: c_upper(kgridN),avg_earnings1
DOUBLE PRECISION   :: avg_earnings_ss(typesN,(shocksN)/2),avg_earnings_ss1(typesN,(shocksN)/2),&
    avg_earnings_ss_ig(typesN,(shocksN)/2)
DOUBLE PRECISION   :: gs,increment_real,increment_real2
DOUBLE PRECISION   :: ss_count(typesN,(shocksN)/2),skink(skinkN),sspay(skinkN)

DOUBLE PRECISION   :: hbar,partprod, K_ig,E,Ec,Ec1,Ep, Ec_ig,N_ig, tr_ig, K, N, N1,N2,K1,K2,Ae,Atfp,sigmaa,gov_spend
DOUBLE PRECISION   :: K1_prime,N1_prime,Y, w, r,kgridl, kgridu,  tr1,lambda_ig,sigmav,sigmaeps,rho
DOUBLE PRECISION   :: a,a1,tax_diff,tax_diff_keep,govt_taxes,tol,pop_start
DOUBLE PRECISION   :: shocks_temp((shocksN)/2),stat_dist_temp((shocksN)/2),shock_pi_temp((shocksN)/2,(shocksN)/2)
DOUBLE PRECISION   :: shocks(shocksN),shock_pi(shocksN,shocksN,Jnr-1,typesN),shock_cumpi(shocksN,shocksN,J,typesN)
DOUBLE PRECISION   :: avg_earnings,avg_earnings_ig,prices(9)
DOUBLE PRECISION   :: random(J,simulationsN),buf(4)

DOUBLE PRECISION   :: labor_income_sorted((Jnr-1)*simulationsN,6)
DOUBLE PRECISION   :: capital_sorted((J+1)*simulationsN,6)



DOUBLE PRECISION, allocatable, dimension(:,:) :: zeros, population,sim_labor_alloc,sim_capital_alloc
DOUBLE PRECISION, allocatable, dimension(:,:) :: sim_consumption_alloc,sim_earnings_alloc,sim_energy_alloc
integer, allocatable, dimension(:) :: zeros2
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: v_alloc,kprime_alloc,labor_alloc


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MPI Parameters

root_process = 0
call MPI_INIT(ierr)
call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
num_procs_run=num_procs
num_procs_run2=num_procs
!Increments for simulations
if (simulationsN<(num_procs-1)) then
    num_procs_run2=simulationsN+1
end if
increment_real2=simulationsN/(num_procs_run2-1.0D0)
increment2=floor(increment_real2)
large_inc2=anint((increment_real2-increment2)*(num_procs_run2-1))

CALL ERSET(4,1,0)
CALL ERSET(3,1,0)
CALL ERSET(2,1,0)

!Increments for value function calculation
if (kgridN<(num_procs-1)) then
    num_procs_run=kgridN+1
end if

increment_real=kgridN/(num_procs_run-1.0D0)
increment=floor(increment_real)
large_inc=anint((increment_real-increment)*(num_procs_run-1))


aggregates=0.0D0
household=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Computational Parameters
a=0.15D0
a1=0.15D0
tol=0.0005D0
!allocate(zeros2(max_outer_loops+1))
!zeros2=0
!taxes_clear=zeros2
!deallocate(zeros2)
taxes_clear=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Utility Parameters
sigma1=2.0D0
if(household==0) then
    !Individual
    sigma2=0.5D0
elseif (household==1) then
    !Household
    sigma2=0.55D0
else
    print*,"Set Household or Individual"
end if

beta_l=0.99502D0
beta_h=0.99502D0
!beta_l=1.0005D0
!beta_h=1.001D0
beta_grid(1)=beta_l
!do parm_it=2,betaN
!    beta_grid(parm_it)=beta_grid(parm_it-1)+(beta_h-beta_l)/real(betaN-1.0D0)
!end do

chi_l=73.1D0
!chi_l=76.0D0
chi_h=73.3D0
chi_grid(1)=chi_l
do parm_it=2,chiN
    chi_grid(parm_it)=chi_grid(parm_it-1)+(chi_h-chi_l)/real(chiN-1.0D0)
end do

!gamma
lambda_l=0.99074D0
!lambda_l=0.99D0
lambda_h=0.99076D0
lambda_grid(1)=lambda_l
do parm_it=2,lambdaN
    lambda_grid(parm_it)=lambda_grid(parm_it-1)+(lambda_h-lambda_l)/real(lambdaN-1.0D0)
end do

!Ebar (defined as theta grid in the code)
theta_l=0.001322D0
!theta_l=0.002D0
theta_h=0.001324D0
theta_grid(1)=theta_l
do parm_it=2,thetaN
    theta_grid(parm_it)=theta_grid(parm_it-1)+(theta_h-theta_l)/real(thetaN-1.0D0)
end do


kgridl_l=-0.157D0
kgridl_h=-0.159D0
kgridl_grid(1)=kgridl_l
do parm_it=2,kgridlN
    kgridl_grid(parm_it)=kgridl_grid(parm_it-1)+(kgridl_h-kgridl_l)/real(kgridlN-1.0D0)
end do




parm_loop=1
do parm_it3=1,lambdaN
do parm_it4=1,thetaN
do parm_it=1,betaN
do parm_it2=1,chiN
do parm_it5=1,1
do parm_it6=1,1
do parm_it7=1,kgridlN
    parms(parm_loop,1)=beta_grid(parm_it)
    parms(parm_loop,2)=chi_grid(parm_it2)
    parms(parm_loop,3)=lambda_grid(parm_it3)
    parms(parm_loop,4)=theta_grid(parm_it4)
    parms(parm_loop,7)=kgridl_grid(parm_it7)
    parm_loop=parm_loop+1
end do
end do
end do
end do
end do
end do
end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Productivity Parameters
!Individual
hbar=0.0D0 !Not used in this code
partprod=1.0D0 !Not used in this code

if(household==0) then
epsilon1(:,1) = (/1.0D0, 1.05037955556458D0, 1.10168574952967D0, 1.15381030928342D0, 1.20663608757434D0, 1.26003734304819D0, 1.3138800980071D0, 1.36802257372016D0, 1.42231570281972D0, 1.47660371749329D0, 1.53072481133652D0, 1.584511871879D0, 1.63779327994293D0, 1.6903937711557D0, 1.74213535412359D0, 1.79283827899664D0, 1.84232204942572D0, 1.89040647024333D0, 1.93691272260084D0, 1.9816644577762D0, 2.02448890043734D0, 2.06521795181501D0, 2.10368928301267D0, 2.13974740856428D0, 2.1732447303485D0, 2.20404254208235D0, 2.23201198484835D0, 2.25703494445744D0, 2.27900488191109D0, 2.29782758879663D0, 2.31342186012357D0, 2.32572007787805D0, 2.33466869942757D0, 2.34022864584003D0, 2.34237558617614D0, 2.34110011486113D0, 2.33640782032644D0, 2.32831924422039D0, 2.31686973160448D0, 2.30210917366396D0, 2.28410164555347D0, 2.26292494305636D0, 2.23867002274592D0, 2.21144035128553D0, 2.18135117038086D0, 2.14852868468884D0/)
elseif(household==1) then
epsilon1(:,1) = (/1.0D0, 1.11972993442023D0, 1.23163188204661D0, 1.33628274802356D0, 1.43423225535758D0, 1.52600294491726D0, 1.61209017543323D0, 1.69296212349824D0, 1.76905978356708D0, 1.84079696795664D0, 1.90856030684589D0, 1.97270924827585D0, 2.03357605814965D0, 2.09146582023247D0, 2.14665643615159D0, 2.19939862539634D0, 2.24991592531814D0, 2.29840469113051D0, 2.34503409590901D0, 2.3899461305913D0, 2.4332556039771D0, 2.47505014272822D0, 2.51539019136855D0, 2.55430901228404D0, 2.59181268572273D0, 2.62788010979474D0, 2.66246300047226D0, 2.69548589158955D0, 2.72684613484296D0, 2.75641389979091D0, 2.7840321738539D0, 2.8095167623145D0, 2.83265628831737D0, 2.85321219286924D0, 2.87091873483891D0, 2.88548299095726D0, 2.89658485581726D0, 2.90387704187394D0, 2.90698507944442D0, 2.90550731670788D0, 2.89901491970559D0, 2.88705187234091D0, 2.86913497637924D0, 2.84475385144809D0, 2.81337093503703D0, 2.77442148249772D0/)
else
    print*,"set household"
end if
!Firm
alpha1 = 0.3D0
alpha2 = 1.0D0-0.403D0
theta=0.03D0
delta = 0.079062D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Heterogeneity Parameters
if(household==0) then
    rho=0.958D0
    sigmaa=0.065D0**0.5D0
    sigmav=(0.017D0/(1-rho**2))**0.5D0
    sigmaeps=0.081D0**0.5D0
elseif(household==1) then
    rho=0.972D0
    sigmaa=0.022D0**0.5D0
    sigmav=(0.012D0/(1-rho**2))**0.5D0
    sigmaeps=0.093D0**0.5D0
else
    print*,"Set Household"
end if

if(my_id==1) then
print*,"got here"
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Calculate Heterogeneity
do i =1,typesN
    types(i)=exp(-(sigmaa)+(real(i-1)/real(typesN-1))*2.0D0*sigmaa)
end do

!    types(1:2)=(/exp(-(sigmaa)), exp(sigmaa)/)
! Stochastic parameters
CALL rouwenhorst(rho,sigmav,shock_pi_temp,shocks_temp,stat_dist_temp,(shocksN)/2)


shocks_temp=exp(shocks_temp)
shocks_temp=shocks_temp/(sum(stat_dist_temp*shocks_temp))

shocks(1:int(shocksN/2))=shocks_temp*exp(-sigmaeps)
shocks(int(shocksN/2)+1:int(shocksN))=shocks_temp*exp(sigmaeps)

shocks_start1=int((((shocksN)/2)+1)/2)
shocks_start2=int(shocks_start1+((shocksN)/2))

!E to E
do it1=1,Jnr-1
do it2=1,typesN
shock_pi(1:int((shocksN)/2),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(1:int((shocksN)/2),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
end do
end do

if(my_id==1) then
print*,"got here2",shocks(1:shocksN)
end if


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Population Parameters
growth = 0.011D0
if(household==0) then
pop_temp1 =(/197316.D0, 197141.D0, 196959.D0, 196770.D0, 196580.D0, 196392.D0, 196205.D0, 196019.D0, 195830.D0, 195634.D0, 195429.D0, 195211.D0, 194982.D0, 194739.D0, 194482.D0, 194211.D0, 193924.D0, 193619.D0, 193294.D0, 192945.D0, 192571.D0, 192169.D0, 191736.D0, 191271.D0, 190774.D0, 190243.D0, 189673.D0, 189060.D0, 188402.D0, 187699.D0, 186944.D0, 186133.D0, 185258.D0, 184313.D0, 183290.D0, 182181.D0, 180976.D0, 179665.D0, 178238.D0, 176689.D0, 175009.D0, 173187.D0, 171214.D0, 169064.D0, 166714.D0, 164147.D0, 161343.D0, 158304.D0, 155048.D0, 151604.D0, 147990.D0, 144189.D0, 140180.D0, 135960.D0, 131532.D0, 126888.D0, 122012.D0, 116888.D0, 111506.D0, 105861.D0, 99957.D0, 93806.D0, 87434.D0, 80882.D0, 74204.D0, 67462.D0, 60721.D0, 54053.D0,	47533.D0,	41241.D0, 35259.D0, 29663.D0, 24522.D0, 19890.D0, 15805.D0, 12284.D0, 9331.D0, 6924.D0, 5016.D0, 3550.D0, 2454.D0 /)
do i=1,(J-1)
   survive(i) = pop_temp1(i+1)/pop_temp1(i)
end do
survive(J) = 0.0D0
adult_equivalence=1.0D0
else
survive=(/1.0D0, 0.999443429520098D0, 0.999453048952094D0, 0.999474759951225D0, 0.999507517658172D0, 0.999550180294631D0, 0.99958688938571D0, 0.999616966794992D0, 0.999637259341625D0, 0.999646435563957D0, 0.999642137784891D0, 0.999638705313936D0, 0.999628264835928D0, 0.999611266702816D0, 0.999587623092204D0, 0.99955347402391D0, 0.999517484919676D0, 0.999474904683543D0, 0.999424216406889D0, 0.999368982231295D0, 0.999311637585569D0, 0.999247262201619D0, 0.99917475362204D0, 0.999096164795857D0, 0.999015366949988D0, 0.99892772761533D0, 0.998832359878012D0, 0.998729478390282D0, 0.998631176268071D0, 0.998535878125485D0, 0.998443264975948D0, 0.998338654631507D0, 0.998221906285934D0, 0.998078357864358D0, 0.997913954710353D0, 0.997725564043014D0, 0.997510605728885D0, 0.997277513805754D0, 0.99702940364917D0, 0.996758720419625D0, 0.996463887451456D0, 0.996136496024502D0, 0.995758248748826D0, 0.995343091935131D0, 0.994874730958277D0, 0.994354016314192D0, 0.993755337681475D0, 0.993065683870379D0, 0.992293047744931D0, 0.991433174915763D0, 0.990457448940283D0, 0.989303075213917D0, 0.987964810901374D0, 0.986475935412401D0, 0.984833077329247D0, 0.982973697702629D0, 0.980769120672249D0, 0.978166137633033D0, 0.975176536129091D0, 0.971754588528722D0, 0.967817497628447D0, 0.963253722024177D0, 0.957955278732916D0, 0.951933505633426D0, 0.945105096525417D0, 0.937431585348977D0, 0.92887640431547D0, 0.919444350437676D0, 0.909124748145856D0, 0.89794551325945D0, 0.885974030114678D0, 0.873304919941224D0, 0.859956433479889D0, 0.842665940409714D0, 0.824214308977887D0, 0.804792055183895D0, 0.785116734289799D0, 0.765756410717599D0, 0.747202819257837D0, 0.729458167443388D0, 0.713669733378468D0/)
survive(J) = 0.0D0
adult_equivalence=(/1.21998398833892D0, 1.26899136580154D0, 1.32056246723579D0, 1.3734138700101D0, 1.42641978111077D0, 1.47860215039009D0, 1.52912105574409D0, 1.57726536021961D0, 1.62244364105098D0, 1.66417539062616D0, 1.70208248938226D0, 1.73588095063074D0, 1.76537293731188D0, 1.79043905067891D0, 1.81103089091152D0, 1.82716388965886D0, 1.83891041451209D0, 1.84639314540638D0, 1.84977872295234D0, 1.84927166869702D0, 1.84510857731432D0, 1.83755258072496D0, 1.8268880841459D0, 1.81341577406915D0, 1.79744789817028D0, 1.77930381714618D0, 1.75930582848249D0, 1.73777526215036D0, 1.71502884823281D0, 1.69137535648057D0, 1.66711250779728D0, 1.64252415765435D0, 1.61787775143512D0, 1.5934220517087D0, 1.5693851374331D0, 1.54597267508807D0, 1.52336646173711D0, 1.50172324001931D0, 1.48117378507042D0, 1.46182226337349D0, 1.44374586353911D0, 1.42699469901482D0, 1.41159198272448D0, 1.39753447363654D0, 1.3847931952624D0, 1.37331442608375D0, 1.36302096190969D0, 1.35381365016324D0, 1.34557319609731D0, 1.33816224094031D0, 1.33142771197095D0, 1.32520344452275D0, 1.31931307591802D0, 1.31357321133145D0, 1.30779686158245D0, 1.30179715285813D0, 1.29539130836514D0, 1.28840490191053D0, 1.28067638341311D0, 1.27206187634362D0, 1.26244024709437D0, 1.25171844627878D0, 1.23983712196005D0, 1.22677650480938D0, 1.21256256519385D0, 1.19727344219358D0, 1.18104614454827D0, 1.16408352353373D0, 1.14666151776743D0, 1.12913666994348D0, 1.11195391549768D0, 1.09896503226485D0, 1.08988645284858D0, 1.0816144615019D0, 1.07430423674565D0, 1.0681188559846D0, 1.06322914430733D0, 1.0598135018135D0, 1.0580577094685D0, 1.05815471348597D0, 1.06030438823754D0/)
end if

dying = 1.D0 - survive
allocate(population(J,2))
population=0.0D0
population(1,1) = 1.0D0
do i=2,(J)
   population(i,1) = survive(i-1)*population((i-1),1)/(1+growth)
end do
pop_start=population(1,1)/sum(population(1:J,1))
population(:,2) = population(:,1)/sum(population(:,1))


if(my_id==1) then
print*,"got here3"
end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Grids


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Social Security Payment Information
b=0.45D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!! Government Parameters
gs=.157D0
govt_taxes=0.0D0
!lambda0=.258D0
lambda1=.031D0
lambdak=.36D0


!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Initial Guesses
lambda_ig=0.836244401D0
K_ig=2.73731287421892D0
N_ig=0.817175505458663D0
Ec_ig=0.01319101D0
!E_ig=16.5373251118981D0
tr_ig=0.0418395847126317D0
ss_ig(2,:)=(/0.247950647D0,0.27492346D0,0.3151075D0,0.341967705D0,0.374516481D0/)
ss_ig(1,:)=(/0.20630098D0,0.210368123D0,0.240250737D0,0.260632005D0,0.289537686D0/)
avg_earnings_ig=0.87436589019486D0
avg_earnings_ss_ig=0.644971442206701D0
avg_earnings_ss_ig(1,:)=(/0.380574867D0,0.393284689D0,0.486667859D0,0.550359321D0,0.640689576D0/)
avg_earnings_ss_ig(2,:)=(/0.510730077D0,0.595020119D0,0.720595242D0,0.804533383D0,0.918847904D0/)

!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_loop = 1
tax_diff=1
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0
if(my_id==1) then
print*,"got here4"
end if

!Set Up population for simulations
seed=53415691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do

i=1
do s=1,typesN
    sim_types(1,i:i+int(simulationsN/(typesN))-1)=s
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(typesN))-1
end if
    i=i+int(simulationsN/(typesN))
end do

!sim_types(1,1:int((simulationsN)/2))=1
!sim_types(1,int((simulationsN)/2)+1:simulationsN)=2
do s=2,J
    sim_types(s,:)=sim_types(1,:)
end do



seed=52345691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do
i=1
do it = 1,typesN
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start1
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start2
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
end do


out_loop=1

!MASTER
if (my_id .eq. root_process) then


do parm_loop=1,betaN*chiN*lambdaN*thetaN*1*1*kgridlN
out_loop=1
beta=parms(parm_loop,1)
chi=parms(parm_loop,2)
gamma1=parms(parm_loop,3)
ebar=parms(parm_loop,4)

kgridl=parms(parm_loop,7)



kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.25D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive

!print*,"kgrid",kgrid
!
!
!
!kgrid(1)=kgridl
!kgrid(kgridnegN)=.0D0
!kgridu = 201.0D0
!kstart = 1
!!Curved grids to handle outliers
!curv = 2.25D0
!
!do ac=2,kgridN
!	kgrid(ac)=kgrid(1)+(kgridu-kgridl)*((ac-1.0D0)/(kgridN-1.0D0))**curv
!end do
!kgrid(kgridN)=kgridu



print*,"in here"

do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
print*,"in here2"

ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it).le.shocksN/2) then
        ss_count(sim_types(1,it),sim_shocks(Jnr,it))=ss_count(sim_types(1,it),sim_shocks(Jnr,it))+1.0D0
    else
        ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+1.0D0
    end if
end do
if(my_id==1) then
    print*,"ss count",ss_count(1,:)
    print*,"ss count",ss_count(2,:)
end if



diff_parameter=1.0D0
Atfp=1.0D0
Ae=Atfp*1.0D0
bad_combination_inner_stop=0
bad_combination_outer=0

!allocate(zeros2(max_outer_loops+1))
!zeros2=0
taxes_clear=0
!deallocate(zeros2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENTER FUNCTION HERE IF MAKING A FUNCTION
!Initial guesses
K=K_ig
N=N_ig
tr=tr_ig
Ec=Ec_ig
!E=E_ig
ss=ss_ig
lambda0=lambda_ig
avg_earnings=avg_earnings_ig
avg_earnings_ss=avg_earnings_ss_ig
!Loop for price convergence
DO WHILE (diff_parameter>tol)
   tax_loop = 1
   tax_diff=1.0D0
   non_clear=0


    passed_parameters(1)=K
    passed_parameters(2)=N
    passed_parameters(3)=Ec
    passed_parameters(4)=Atfp
    passed_parameters(5)=Ae
    passed_parameters(6)=0.0D0

    if(out_loop==1) then
        xguess(1)=0.05D0
    else
        xguess(1)=Ep
    end if
    if(out_loop==1) then
        xguess(2)=K*.5D0
        xguess(3)=N*.6D0
    else
        xguess(2)=K*K1_share
        xguess(3)=N*N1_share
    end if

    print*,"parms",Atfp,Ae,xguess,Ec
    print*,"guess",xguess
    print*,"parts",N,K
    CALL D_NEQNF(energy_root_finder,solver_out,ERRREL=0.0001D0,xguess=xguess)
    Ep=solver_out(1)
    E=Ep+Ec
    K1=solver_out(2)
    K2=K-K1
    N1=solver_out(3)
    N2=N-N1
    K1_share=K1/K
    N1_share=N1/N

    Y=Atfp*(K1**(alpha1)*N1**(1.0D0-alpha1-theta))*Ep**theta
    pe=Atfp*theta*(K1/Ep)**alpha1*(N1/Ep)**(1.0D0-alpha1-theta)
    gov_spend=(Y+pe*Ec)*gs
    w=Atfp*(1.0D0-alpha1-theta)*(K1/N1)**(alpha1)*(Ep/N1)**theta
    r=Atfp*alpha1*(N1/K1)**(1.0D0-alpha1-theta)*(Ep/K1)**theta-delta
    tauss=sum(sum(ss*ss_count,2),1)/real(simulationsN)*sum(population(Jnr:J,2))/&
        (sum(sum(avg_earnings_ss*ss_count,2),1)/real(simulationsN)*sum(population(1:Jnr-1,2)))
    ss_cap=avg_earnings*2.42
print*,"pe Y w r",pe,Y,w,r,N,Ep
    avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)
!
!Extra loop for convergence on tax parameters (if turned on need tax parameter convergence begore prices update)
    DO WHILE (tax_diff >tol)
!Solve for policy funcation
        v = -1.0D0/0.0D0
        v(:,:,:,J+1)=0.0D0
        v_retired = -1.0D0/0.0D0
        kprime=0
        v_retired(:,:,capital_positive:kgridN,J+1)=0.0D0

!Last generation
            do it=capital_positive,kgridN
                do it2=1,typesN
                    do e_it=1,shocksN
                        if(e_it .le. (shocksN)/2) then
                            ctot80(it,1)=retcon(kgrid(it),0.0D0,tr,r,ss(it2,e_it),avg_earnings)
                            energy_miss=energy_sharer(0.00D0,ctot80(it,1)*.15D0,ctot80(it,1),energy_share_equation,.0001D0,c80(it,1),ctot80(it,1),J)
                            e80(it,1)=(ctot80(it,1)-c80(it,1))/pe
                            v_retired(it2,e_it,it,J) =utility_last(c80(it,1),e80(it,1),adult_equivalence(J))
                        else
                            v_retired(it2,e_it,it,J) =v_retired(it2,e_it-(shocksN/2),it,J)
                        end if
                    end do
                end do
            end do
!print*,"v_retired",s,v_retired(2,:,20,J),capital_positive
!print*,"v_retired",J,v_retired(1,1,capital_positive-10:capital_positive+30,J)

!Retired generations
        allocate(zeros(kgridN,1))
        zeros=0.0D0
        do s =(J-1),Jnr, -1
            do e_it=1,shocksN
            if(e_it .le. (shocksN)/2) then
                do it =1,typesN
                do kit=capital_positive, kgridN
                    con=-1.0D0
                    con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,kgrid(kit),kgrid(capital_positive:),&
                        tr,r,ss(it,e_it),avg_earnings)
                    con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                    if (con(kgridN)>0.015+ebar) then
                        con_pos=kgridN
                    end if
                    if (con(capital_positive+1)<=0.015D0+ebar) then
                        kprime(it,e_it,kit,s)=0
                        v_retired(it,e_it,kit,s)=-1.0D0/0.0D0
                    else
                    CALL retired_valuer(kgrid(kit),0.0D0,kgrid(con_pos),ss(it,e_it),tr,r,survive(s),v_retired(it,e_it,:,s+1),&
                        v_retired(it,e_it,kit,s),kprime(it,e_it,kit,s),avg_earnings,s)
                    end if
               end do
               end do
            else
                do it=1,typesN
                do kit=capital_positive, kgridN
                    kprime(it,e_it,kit,s)=kprime(it,e_it-(shocksN/2),kit,s)
                    v_retired(it,e_it,kit,s)=v_retired(it,e_it-(shocksN/2),kit,s)
                end do
               end do
            end if
            end do
!print*,"v_retired",s,v_retired(2,:,21,s)
!print*,"v_retired",s,v_retired(1,1,capital_positive-10:capital_positive+30,s)
!print*,"v_retired",s,v_retired(2,:,20,s)
        end do
        deallocate(zeros)
!        do e_it=1,shocksN
!            do i=1,typesN
!                v_retired(i,e_it,:,Jnr:J-1)=v_retired(1,1,:,Jnr:J-1)
!                    kprime(i,e_it,:,Jnr:J-1)=kprime(1,1,:,Jnr:J-1)
!            end do
!        end do



!Working generations
v(:,:,:,Jnr)=v_retired(:,:,:,Jnr)
                prices(1)=tr
                prices(2)=lambda0
                prices(3)=w
                prices(4)=r
                prices(5)=tauss
!                prices(6)=ss(1)
                prices(7)=ss_cap
                prices(8)=avg_earnings
                prices(9)=pe

                CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
                CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
        do s =(Jnr-1),2,-1
            CALL MPI_Bcast( v(:,:,:,s+1),typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            end_row=0
            do an_id = 1, num_procs_run -1
                start_row=end_row+1
                if (an_id>large_inc) then
                    end_row=start_row+increment-1
                else
                    end_row=start_row+increment
                end if
                call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
                call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
            end do
            do an_id=1,num_procs_run -1
                call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( labor(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( kprime(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( v(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            end do
!print*,"v",s,v(2,:,40,s)
!print*,"kprime",s,kprime(2,:,40,s)
!print*,"labor",s,labor(2,:,40,s)
        end do




!! Simulate peoples
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end_row=0
    do an_id = 1, num_procs_run2 -1
        start_row=end_row+1
        if (an_id>large_inc2) then
            end_row=start_row+increment2-1
        else
            end_row=start_row+increment2
        end if
        call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
        call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
    end do
    do an_id=1,num_procs_run2 -1
        call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_labor(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_earnings(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_capital(:,start_row:end_row),(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_consumption(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_energy(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    end do








!Calculate Aggregate Profiles

agg_labor(:,1)=sum(sim_labor(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_labor(:,2)=sum(sim_labor(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_capital(:,1)=sum(sim_capital(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_capital(:,2)=sum(sim_capital(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_consumption(:,1)=sum(sim_consumption(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_consumption(:,2)=sum(sim_consumption(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_energy(:,1)=sum(sim_energy(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_energy(:,2)=sum(sim_energy(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
!
!print*,"sim consumption1",agg_consumption(:,1)
!print*,"sim consumption2",agg_consumption(:,2)
!
!print*,"sim labor1",agg_labor(:,1)
!print*,"sim labor2",agg_labor(:,2)
!
!print*,"sim savings1",agg_capital(:,1)
!print*,"sim savings2",agg_capital(:,2)
!print*,"sim labor",sim_labor(:,1)

!!Find lambda 2 that clears the market
CALL lambda2_optimizer(lambda0_1,gov_spend,sim_labor,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)



        if (lambda0_1>15)then
            a1=0.20D0
        else
            a1=0.20D0
        end if
        tax_diff=abs((lambda0_1)-lambda0)/lambda0
        tax_diff_keep=tax_diff
        lambda0=(a1*lambda0_1+(1-a1)*lambda0)

        govt_taxes=tax_differ2(lambda0,sim_labor,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend,population(:,2),avg_earnings,ss,survive)
        !Check for bad convergence within tax loop
        if (value==-1000000) then
            tax_diff=0
            bad_combination_inner_stop=1
        end if

        tax_loop=tax_loop+1
        if (tax_loop > max_tax_loop) then
          tax_diff = 0
          bad_combination_inner_nostop=1
        end if
        CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    END DO




    K1_prime=0.0D0
    N1_prime=0.0D0
    tr1=0.0D0
    hours_jnr=0.0D0
    Ec1=0.0D0
    N1_ss=0.0D0

    avg_earnings1=0.0D0
    avg_earnings_ss1=0.0D0
    borrowers=0.0D0
    negative_value=0.0D0
    negative_count=0.0D0
    !Calculate new aggregates
    do it=1,simulationsN
        do s=1,Jnr-1
            avg_earnings1=avg_earnings1+sim_earnings(s,it)*population(s,2)/(simulationsN*sum(population(1:Jnr-1,2)))
            if (sim_shocks(Jnr,it) .le. shocksN/2) then
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it))*&
                    sum(population(1:Jnr-1,2)))
            else
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)*&
                    sum(population(1:Jnr-1,2)))
            end if
        end do
        do s=1,Jnr-1
            if (sim_labor(s,it)<hbar) then
                N1_prime=N1_prime+partprod*sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
                N1_ss(sim_types(s,it))=N1_ss(sim_types(s,it))+partprod*sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN*real(typesN)
            else
                N1_prime=N1_prime+sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
                N1_ss(sim_types(s,it))=N1_ss(sim_types(s,it))+sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN*real(typesN)
            end if
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
            hours_jnr=hours_jnr+sim_labor(s,it)*population(s,2)/simulationsN
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
            if(sim_capital(s,it)<0.0D0) then
                borrowers=borrowers+population(s,2)/simulationsN
                negative_count=negative_count+population(s,2)/simulationsN
                negative_value=negative_value+sim_capital(s,it)*population(s,2)/simulationsN
            end if

        end do
        do s=Jnr,J
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
            if(sim_capital(s,it)<0.0D0) then
                borrowers=borrowers+population(s,2)/simulationsN
                negative_count=negative_count+population(s,2)/simulationsN
                negative_value=negative_value+sim_capital(s,it)*population(s,2)/simulationsN
            end if
        end do
    end do

        negative_value=negative_value/negative_count/(K)
        hours_jnr=hours_jnr/sum(population(1:Jnr-1,2))
        K1_prime=K1_prime/(1.0D0+growth)
        tr1=tr1/(1.0D0+growth)

    diff_parameter_v(1)=abs(((tr1)-tr))/tr
    diff_parameter_v(2)=abs(((K1_prime)-K))/K
    diff_parameter_v(3)=abs(((N1_prime)-N))/N
!    diff_parameter_v(4)=abs(((ss1*a+ss*(1-a))-ss))/ss
!    ,abs(((ss1(2))-ss(2)))/ss(2))
    diff_parameter_v(5)=abs(tax_diff_keep)
    diff_parameter_v(6)=max(abs(((avg_earnings1)-avg_earnings))/avg_earnings,abs(((Ec1)-Ec))/Ec,&
        maxval(abs(((avg_earnings_ss1)-avg_earnings_ss))/avg_earnings_ss1))
!    diff_parameter=MAXVAL(diff_parameter_v)
    print*,"K1",K1_prime,K
    print*,"N1",N1_prime,N
    print*,"tr1",tr1,tr
    print*,"hours",hours_jnr
    print*,"avg earnings1",avg_earnings_ss1(1,:)
    print*,"avg earnings2",avg_earnings_ss1(2,:)

!!!!Update Aggregates
    Ec=(Ec1*a+Ec*(1-a))
    N=(N1_prime*a+N*(1-a))
    K=(K1_prime*a+K*(1-a))
!    do it=1,typesN
!    end do
    tr=(tr1*a+tr*(1.0D0-a))
    avg_earnings=(avg_earnings1*a+avg_earnings*(1.0D0-a))
    avg_earnings_ss=(avg_earnings_ss1*a+avg_earnings_ss*(1.0D0-a))



    skink(1)=0.0D0
    skink(2)=0.21D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(3)=1.29D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(4)=2.42D0*sum(avg_earnings_ss*ss_count)/sum(ss_count)
    skink(5)=10000000.0D0*1.4D0
!    if (sgridu*1.4D0<2.42D0*avg_earnings) then
!    skink(5)=avg_earnings*3.0D0
!    end if
    sspay(1)=0.0D0
    sspay(2)=skink(2)*0.90D0
    sspay(3)=sspay(2)+(skink(3)-skink(2))*0.32D0
    sspay(4)=sspay(3)+(skink(4)-skink(3))*0.15D0
    sspay(5)=sspay(4)

print*,"skink sspay",skink,sspay

!    do it=1,sgridN
!        spayoutgrid(it)=spayer(skinkN,skink,sspay,sgrid(it))*retire_scaler(retire_age)
!    end do
    do e_it=1,shocksN/2
        do it=1,typesN
            ss1(it,e_it)=spayer(skinkN,skink,sspay,avg_earnings_ss1(it,e_it))
        end do
    end do
    diff_parameter_v(4)=maxval(abs(((ss1)-ss))/ss)
    print*,"diff_parmaeter",diff_parameter_v
    diff_parameter=MAXVAL(diff_parameter_v)
    print*,"ss1",ss1(1,:),ss1(2,:)

    do e_it=1,shocksN/2
        do it=1,typesN
            ss(it,e_it)=(ss1(it,e_it)*a+ss(it,e_it)*(1-a))
        end do
    end do


print*,"got here 2"

    if (bad_combination_inner_stop == 1) then
        diff_parameter=0.0D0
    end if
print*,"got here 2aa"
    taxes_clear(out_loop)=0
print*,"got here 2a"
    if (abs(govt_taxes/gov_spend)>tol) then
!        taxes_clear(out_loop)=1
        non_clear=1
        diff_parameter=1.0D0
    end if
print*,"got here 2b"
    if (taxes_clear(out_loop)==0) then
        allocate(zeros2(max_outer_loops+1))
        zeros2=0
        taxes_clear=zeros2
        deallocate(zeros2)
    end if
print*,"got here 2"
    if (sum(taxes_clear)>max_non_clear) then
        diff_parameter=0.0D0
        bad_combination_inner_stop=1
    end if
print*,"got here 2d"
    if (out_loop > max_outer_loops) then
        bad_combination_outer=1
        diff_parameter=0.0D0
    end if

    energy_share_sim=0.0D0
    do it=1,simulationsN
        do s=1,J
            energy_share_sim(s,it)=sim_energy(s,it)*pe/((sim_energy(s,it)*pe)+sim_consumption(s,it))
        end do
    end do

    out_loop=out_loop+1
    print*,"Loop"
    print*,"Plus one"
    print*,out_loop,parm_loop
    print*,"lambda0",lambda0
    print*,"govt miss",govt_taxes/gov_spend
    print*,"diffparm",diff_parameter
    print*,"taxdiff",tax_diff_keep
    print*,"hours",hours_jnr
    print*,"K",K
    print*,"N",N
    print*,"Y",Y
    print*,"ss",ss
    print*,"tr",tr

    print*,"max capital",maxval(sim_capital)
    print*,bad_combination_inner_nostop
    CALL MPI_Bcast( diff_parameter, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )

print*,"got here"

do s=1,J
    consumption_temp_sorted=sim_consumption(s,:)+sim_energy(s,:)*pe
    CALL QsortC(consumption_temp_sorted)

    e_c_l=1
    e_c_h=int((simulationsN)/2)
    do it=1,simulationsN
        if(sim_consumption(s,it)+sim_energy(s,it)*pe.le.consumption_temp_sorted(int((simulationsN)/2)) .and. e_c_l<int((simulationsN)/2))then
            energy_temp(e_c_l)=energy_share_sim(s,it)
            e_c_l=e_c_l+1
        else
            energy_temp(e_c_h)=energy_share_sim(s,it)
            e_c_h=e_c_h+1
        end if
    end do
    energy_high(s)=sum(energy_temp(int((simulationsN)/2)+1:simulationsN),1)/real(int((simulationsN)/2))
    energy_low(s)=sum(energy_temp(1:int((simulationsN)/2)),1)/real(int((simulationsN)/2))
end do

print*,"got here1"

END DO

print*,"got here2"

do s=1,J
    consumption_temp_sorted=sim_consumption(s,:)+sim_energy(s,:)*pe
    CALL QsortC(consumption_temp_sorted)

    e_c_l=1
    e_c_h=int((simulationsN)/2)
    do it=1,simulationsN
        if(sim_consumption(s,it)+sim_energy(s,it)*pe.le.consumption_temp_sorted(int((simulationsN)/2)) .and. e_c_l<int((simulationsN)/2))then
            energy_temp(e_c_l)=energy_share_sim(s,it)
            consumption_temp(e_c_l)=sim_consumption(s,it)+sim_energy(s,it)*pe
            e_c_l=e_c_l+1
        else
            energy_temp(e_c_h)=energy_share_sim(s,it)
            consumption_temp(e_c_h)=sim_consumption(s,it)+sim_energy(s,it)*pe
            e_c_h=e_c_h+1
        end if
    end do
    energy_high(s)=sum(energy_temp(int((simulationsN)/2)+1:simulationsN),1)/real(int((simulationsN)/2))
    energy_low(s)=sum(energy_temp(1:int((simulationsN)/2)),1)/real(int((simulationsN)/2))
    consumption_high(s)=sum(consumption_temp(int((simulationsN)/2)+1:simulationsN),1)/real(int((simulationsN)/2))
    consumption_low(s)=sum(consumption_temp(1:int((simulationsN)/2)),1)/real(int((simulationsN)/2))
end do


temp=1
do i=1,simulationsN
    labor_income_sorted(temp:temp+(Jnr-1)-1,1)=sim_earnings(1:Jnr-1,i)
    labor_income_sorted(temp:temp+(Jnr-1)-1,2)=population(1:Jnr-1,2)
    temp=temp+Jnr-1
end do

do it= 1, (Jnr-1)*simulationsN
    krow = minloc( labor_income_sorted( it:(Jnr-1)*simulationsN, 1 ), dim=1 ) + it - 1
    buf( 1: 2)     = labor_income_sorted( it, 1:2 )
    labor_income_sorted( it, 1:2 ) = labor_income_sorted( krow, 1:2 )
    labor_income_sorted( krow, 1:2 ) = buf( 1:2 )
    if(it>1) then
        labor_income_sorted(it,3)=labor_income_sorted(it-1,3)+labor_income_sorted(it,1)*labor_income_sorted(it,2)
        labor_income_sorted(it,4)=labor_income_sorted(it-1,4)+labor_income_sorted(it,2)
    else
        labor_income_sorted(it,3)=labor_income_sorted(it,1)*labor_income_sorted(it,2)
        labor_income_sorted(it,4)=labor_income_sorted(it,2)
    end if
end do
labor_income_sorted(:,3)=labor_income_sorted(:,3)/labor_income_sorted(simulationsN*(Jnr-1),3)
labor_income_sorted(:,4)=labor_income_sorted(:,4)/labor_income_sorted(simulationsN*(Jnr-1),4)

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.99))
print*,"labor income 1",labor_income_sorted(intloc1,1),1.0D0-labor_income_sorted(intloc1,3),intloc1

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.95))
print*,"labor income 5",labor_income_sorted(intloc1,1),1.0D0-labor_income_sorted(intloc1,3),intloc1

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.80))
print*,"labor income 20",labor_income_sorted(intloc1,1),1.0D0-labor_income_sorted(intloc1,3),intloc1

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.60))
print*,"labor income 40",labor_income_sorted(intloc1,1),1.0D0-labor_income_sorted(intloc1,3),intloc1

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.20))
print*,"labor income bottom 20",labor_income_sorted(intloc1,1),labor_income_sorted(intloc1,3),intloc1

intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.10))
print*,"labor income bottom 10",labor_income_sorted(intloc1,1),labor_income_sorted(intloc1,3),intloc1

temp=1
do i=1,simulationsN
    capital_sorted(temp:temp+(J+1)-1,1)=sim_capital(1:J+1,i)
    capital_sorted(temp:temp+(J+1)-1,2)=population(1:J+1,2)
    temp=temp+J+1
end do

do it= 1, (J+1)*simulationsN
    krow = minloc( capital_sorted( it:(J+1)*simulationsN, 1 ), dim=1 ) + it - 1
    buf( 1: 2)     = capital_sorted( it, 1:2 )
    capital_sorted( it, 1:2 ) = capital_sorted( krow, 1:2 )
    capital_sorted( krow, 1:2 ) = buf( 1:2 )
    if(it>1) then
        capital_sorted(it,3)=capital_sorted(it-1,3)+capital_sorted(it,1)*capital_sorted(it,2)
        capital_sorted(it,4)=capital_sorted(it-1,4)+capital_sorted(it,2)
    else
        capital_sorted(it,3)=capital_sorted(it,1)*capital_sorted(it,2)
        capital_sorted(it,4)=capital_sorted(it,2)
    end if
end do
capital_sorted(:,3)=capital_sorted(:,3)/capital_sorted(simulationsN*(J+1),3)
capital_sorted(:,4)=capital_sorted(:,4)/capital_sorted(simulationsN*(J+1),4)

!print*,"got here3"


aggregates(parm_loop,1)=K/(Y+pe*Ec)
aggregates(parm_loop,2)=hours_jnr
aggregates(parm_loop,3)=K
aggregates(parm_loop,4)=N
aggregates(parm_loop,5)=tr
aggregates(parm_loop,6)=lambda0
aggregates(parm_loop,7)=gov_spend
aggregates(parm_loop,8)=tauss
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.80))
aggregates(parm_loop,9)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.80))
aggregates(parm_loop,10)=1.0D0-capital_sorted(intloc1,3)
aggregates(parm_loop,11)=sum(sum(energy_share_sim,2)*population(1:J,2),1)/real(simulationsN)
aggregates(parm_loop,12)=sum(energy_low*population(1:J,2))
aggregates(parm_loop,13)=sum(energy_high*population(1:J,2))
aggregates(parm_loop,14)=avg_earnings
aggregates(parm_loop,15)=E
aggregates(parm_loop,16)=beta
aggregates(parm_loop,17)=chi
aggregates(parm_loop,18)=gamma1
aggregates(parm_loop,19)=ebar
aggregates(parm_loop,20)=kgridl
aggregates(parm_loop,21)=K1
aggregates(parm_loop,22)=N1
aggregates(parm_loop,23)=out_loop
aggregates(parm_loop,24)=0.0D0
aggregates(parm_loop,25)=Y
aggregates(parm_loop,26)=pe
aggregates(parm_loop,27)=Ep
aggregates(parm_loop,28)=0.0D0
aggregates(parm_loop,29)=sum(consumption_low*population(1:J,2))
aggregates(parm_loop,30)=sum(consumption_high*population(1:J,2))
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.99))
aggregates(parm_loop,31)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.95))
aggregates(parm_loop,32)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.90))
aggregates(parm_loop,33)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.80))
aggregates(parm_loop,34)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.60))
aggregates(parm_loop,35)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.40))
aggregates(parm_loop,36)=1.0D0-labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.20))
aggregates(parm_loop,37)=labor_income_sorted(intloc1,3)
intloc1=int(MINLOC(labor_income_sorted(:,4),1,mask=labor_income_sorted(:,4).gt.0.10))
aggregates(parm_loop,38)=labor_income_sorted(intloc1,3)

intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.99))
aggregates(parm_loop,39)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.95))
aggregates(parm_loop,40)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.90))
aggregates(parm_loop,41)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.80))
aggregates(parm_loop,42)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.60))
aggregates(parm_loop,43)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.40))
aggregates(parm_loop,44)=1.0D0-capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.20))
aggregates(parm_loop,45)=capital_sorted(intloc1,3)
intloc1=int(MINLOC(capital_sorted(:,4),1,mask=capital_sorted(:,4).gt.0.10))
aggregates(parm_loop,46)=capital_sorted(intloc1,3)
aggregates(parm_loop,47)=borrowers
aggregates(parm_loop,48)=negative_value

print*,"got here4"

aggregates_ss(parm_loop,:,:)=ss
print*,"got here5"

aggregates_ss_earnings(parm_loop,:,:)=avg_earnings_ss


print*,"out",aggregates(parm_loop,:)
print*,"out1",ss,aggregates_ss(parm_loop,:,:)
print*,"out2",avg_earnings_ss,aggregates_ss_earnings(parm_loop,:,:)
Ec_ig=Ec
K_ig=K
N_ig=N
tr_ig=tr
ss_ig=ss
lambda_ig=lambda0
avg_earnings_ig=avg_earnings
avg_earnings_ss_ig=avg_earnings_ss
end do
if(my_id==1 .and. parm_loop==1) then
print*,"labor",agg_labor(:,1)
print*,"capital",agg_capital(:,1)
print*,"consumption",agg_consumption(:,1)
end if
 OPEN(UNIT=19, FILE="/<filepath calibration>/aggregates_output.dat",STATUS="replace", &
       FORM="formatted")
    do it2=1,betaN*chiN*thetaN*lambdaN*1*1*kgridlN
  WRITE(19,887) aggregates(it2,:)
  end do
  CLOSE(UNIT=13)
887 format(48F24.16)
886 format(i9)
889 format(2F24.16)
1887 format(F24.16)

 OPEN(UNIT=19, FILE="/<filepath calibration>/aggregates_ss.dat",STATUS="replace", &
       FORM="formatted")
    do it2=1,betaN*chiN*thetaN*lambdaN*1*1*kgridlN
  WRITE(19,1887) aggregates_ss(it2,:,:)
  end do
  CLOSE(UNIT=13)

 OPEN(UNIT=19, FILE="/<filepath calibration>/ss_earnings_output.dat",STATUS="replace", &
       FORM="formatted")
    do it2=1,betaN*chiN*thetaN*lambdaN*1*1*kgridlN
  WRITE(19,1887) aggregates_ss_earnings(it2,:,:)
  end do
  CLOSE(UNIT=13)



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!Slaves that are going to work
else if (my_id < (num_procs_run)) then

do parm_loop=1,betaN*chiN*lambdaN*thetaN*1*1*kgridlN
out_loop=1
beta=parms(parm_loop,1)
chi=parms(parm_loop,2)
gamma1=parms(parm_loop,3)
ebar=parms(parm_loop,4)
kgridl=parms(parm_loop,7)


kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.25D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive






do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do

ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it)>shocksN/2) then
        ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+1.0D0
    else
        ss_count(sim_types(1,it),sim_shocks(Jnr,it))=ss_count(sim_types(1,it),sim_shocks(Jnr,it))+1.0D0
    end if
end do
if(my_id==1) then
    print*,"ss count",ss_count(1,:)
    print*,"ss count",ss_count(2,:)
end if


diff_parameter=1.0D0
Atfp=1.0D0
Ae=Atfp*1.0D0

DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)
            CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
!            ss=prices(6)
            ss_cap=prices(7)
            avg_earnings=prices(8)
            pe=prices(9)
            avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)



!Solve Working
        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            call MPI_RECV( end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            allocate(labor_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(kprime_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(v_alloc(typesN,shocksN,end_row-start_row+1))
            labor_alloc=0.0D0
            kprime_alloc=0
            v_alloc=0.0D0
            do kit=start_row,end_row
                if(kgrid(kit)<0.0D0) then
                    rate_use=r/survive(s-1)
                else
                    rate_use=r
                end if
                do e_it=1,shocksN
                    do i=1,typesN
                        c_upper=c_upper_fnc(epsilon1(s,1),kgrid(kit),kgrid,kgridN,tr,&
                            w*shocks(e_it)*types(i),rate_use)
                        temp=int(MINLOC(c_upper,1,mask=c_upper.gt.0.015D0+ebar))
                        if (c_upper(kgridN)>0.015D0+ebar) then
                            temp = kgridN
                        end if
                        if (c_upper(2)<=0.015D0+ebar) then
                            kprime_alloc(i,e_it,kit-start_row+1)=0
                            v_alloc(i,e_it,kit-start_row+1)=-1.0D0/0.0D0
                        else
                            CALL working_valuer(s,kgrid(kit),temp,tr,rate_use,&
                            survive(s),&
                            w*shocks(e_it)*types(i)*epsilon1(s,1),shock_pi(e_it,:,s,i),&
                            v_in(i,:,1:temp),&
                            v_alloc(i,e_it,kit-start_row+1),kprime_alloc(i,e_it,kit-start_row+1),&
                            labor_alloc(i,e_it,kit-start_row+1),&
                            partprod,hbar)

!if(i==1 .and. e_it==1 .and. kit==40) then
!print*,"working guts",v_alloc(i,e_it,kit-start_row+1),kprime_alloc(i,e_it,kit-start_row+1),&
!                            labor_alloc(i,e_it,kit-start_row+1)
!print*,"value prime",temp,v_in(i,1,35:45)
!
!end if
                        end if
                    end do
                end do
            end do
            call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( labor_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( kprime_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( v_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            deallocate(labor_alloc)
            deallocate(kprime_alloc)
            deallocate(v_alloc)
        end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))



    !!Initialize
    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0





    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
        end do


        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if
            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do
end do



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves that are not going to work for some
else if (my_id > (num_procs_run-1)) then
do parm_loop=1,betaN*chiN*lambdaN*thetaN*1*1*kgridlN
out_loop=1
beta=parms(parm_loop,1)
chi=parms(parm_loop,2)
gamma1=parms(parm_loop,3)
ebar=parms(parm_loop,4)
kgridl=parms(parm_loop,7)


kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.25D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive





do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do

ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it)>shocksN/2) then
        ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+1.0D0
    else
        ss_count(sim_types(1,it),sim_shocks(Jnr,it))=ss_count(sim_types(1,it),sim_shocks(Jnr,it))+1.0D0
    end if
end do
if(my_id==1) then
    print*,"ss count",ss_count(1,:)
    print*,"ss count",ss_count(2,:)
end if

diff_parameter=1.0D0
Atfp=1.0D0
Ae=Atfp*1.0D0

DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving Value Function (These are not working)

            CALL MPI_Bcast( prices, 9, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
!            ss=prices(6)
            ss_cap=prices(7)
            avg_earnings=prices(8)
            pe=prices(9)
            avg_earnings_scale=avg_earnings*(1.0D0-tauss/2.0D0)



        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
       end do




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))



    !!Initialize

    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0





    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
        end do


        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it)-shocksN/2,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if
            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/pe
        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do

end do




end if
deallocate(population)








call MPI_FINALIZE(ierr)

end PROGRAM hetero_ss
